The following prior art references, which are incorporated by reference in their entirety herein, are considered to be relevant for an understanding of the invention:                [1] A. Bamberger, P. Joly and J. E. Roberts. Second-order absorbing boundary conditions for the wave equation: A solution to the corner problem. SIAM J. on Numerical Analysis 27(2):323-352, 1990.        [2] A. Bayliss, M. Gunzburger and E. Turkel. Boundary Conditions for the numerical solution of elliptic equations in exterior regions. SIAM J. of Applied Mathematics 42:430-451, 1982.        [3] J.-P. Berenger. A perfectly matched layer for the absorption of electromagnetic waves. J. of Computational Physics 114:185-200, 1994.        [4] J.-P. Berenger. Three dimensional perfectly matched layer for the absorption of electromagnetic waves. J. of Computational Physics 127:363-379, 1996.        [5] B. Engquist and A. Majda. Absorbing boundary conditions for the numerical simulation of waves. Mathematics of Computation 31:629-651, 1977.        [6] D. Givoli. High-order local non-reflecting boundary conditions: a review. Wave Motion 39(4):319-326, 2004.        [7] D. Gordon and R. Gordon. CARP-CG: a robust and efficient parallel solver for linear systems, applied to strongly convection dominated PDEs. Parallel Computing 36:495-515, 2010.        [8] Y. Li, L. Metivier, R. Brossier, B. Han and J. Virieux. 2D and 3D frequency-domain elastic wave modeling in complex media with a parallel iterative solver. Geophysics 80(3):T101T118, 2015.        [9] M. Medvinsky, E. Turkel and U. Hetmaniuk. Local absorbing boundary conditions for elliptical shaped bodies. J. of Computational Physics 227:8254-8267, 2008.        [10] J. Pujol. Elastic Wave Propagation and Generation in Seismology. Cambridge University Press, Cambridge, UK, 2003.        [11] A. Sommerfeld. Partial Differential Equations in Physics, Section 28, pages 188-200, Academic Press, New York, 1964.        [12] E. Turkel, D. Gordon, R. Gordon and S. Tsynkov. Compact 2D and 3D sixth order schemes for the Helmholtz equation with variable wave Number. J. of Computational Physics 232(1):272-287, 2013.        [13] T. van Leeuwen, D. Gordon, R. Gordon, and F. Herrmann. Preconditioning the Helmholtz equation via row projections. In Proc. 74th EAGE Conference, Copenhagen, Denmark, 60-65. EAGE, June 2012.        [14] T. van Leeuwen and F. Herrmann. 3D frequency-domain seismic inversion with controlled sloppiness. SIAM J. on Scientific Computing 36(5):5192-5217, 2014.        [15] J. Virieux and S. Operto. An overview of full-waveform inversion in exploration geophysics. Geophysics 74(6):wcc127wcc152, 2009.        [16] D. Gordon, R. Gordon, and E. Turkel Compact high order schemes with gradient-direction derivatives for absorbing boundary conditions., Journal of Computational Physics, Volume 297, 15 Sep. 2015, Pages 295-315. (Not prior art as it was published after the filing date of this application and is a publication of the present inventors).        [17] U.S. Pat. No. 6,687,659 to Y. Shen.        
Many physical phenomena are governed by partial differential equations (PDEs). Examples of such phenomena are acoustic waves (governed by the wave equation), electromagnetic radiation (Maxwell's equations), earthquakes (the elastic wave equation), and the quantum state of a physical system (Schrödinger's equation). The propagation of such physical phenomena through three-dimensional (3D) space is often caused by some physical occurrence; for example, earthquakes create seismic waves which propagate in a large volume of the earth. For such physical phenomena, there is often a practical need to obtain information about the material in which the phenomenon occurred. A case in point is exploration geophysics, the science lying behind various techniques for obtaining a three-dimensional image of a volume of the earth. Such information is important in order to detect the composition of the volume and pinpoint important areas of special interest such as water or oil. In the case of exploration geophysics, raw data is typically obtained by sensors placed over the surface of interest, and special machinery impacts the surface at many points. The sensors measure the vibrations caused by the impacts, and the data is gathered and prepared for processing by electronic computers. There are several methods for obtaining an image of the volume of interest. One of the modern methods is called “full waveform inversion”, abbreviated FWI. This method is well-known in the field, and it is described in detail in [15]. A crucial step in FWI, which is repeated many times, is the solution of the Helmholtz equation (also known as the “wave equation in the frequency domain”). To this end, the data is modeled by a three-dimensional grid of points corresponding to the volume of interest. The input data to the solver of the Helmholtz equation consists of the value of the speed of sound at every grid point. A numerical method is then applied to the data with the purpose of solving the Helmholtz equation for the given data set. The solution of the Helmholtz equation is a displacement value at every grid point. For the purpose of expediting the calculations, it is essential to perform the calculations as efficiently as possible so as to minimize the number of required computers, their memory, and the time to complete the calculations. In addition, it is essential to obtain as accurate a solution as possible for the Helmholtz equation.
The reverberations in the earth generated by the impacts travel well beyond the actual region of interest. However, the grid imposed by the numerical computations must be finite, because the computer cannot handle a potentially infinite domain. As a result, it is necessary to impose certain conditions on the boundary of the grid that will simulate the fact that the reverberations actually travel beyond the region of interest. Such numerical conditions are called “absorbing boundary conditions” (ABCs), or “non-reflecting boundary conditions”. In addition to the numerical conditions, some ABCs also require certain conditions to be imposed on the physical configuration of the problem being solved. For example, there may be restrictions on the geometrical shape of the domain and/or the position of the source of impact causing the physical phenomenon within the domain.
The theory lying behind ABCs is very extensive and generally well known by practitioners of the art. For a review of ABCs, see [6]. ABCs are used to impose certain numerical conditions on the boundaries of the finite grid which correspond to the physical domain of interest. ABCs simulate the fact that in the real world, physical phenomena caused by some disturbance in the domain of interest have an effect well beyond the relatively small domain. Without ABCs, the obtained numerical solution can appear to be reflected from the boundaries. Accurate ABCs are essential to the process of solving the Helmholtz equation; without them, the solution of the equation may not be sufficiently accurate, and this might hamper the entire FWI process. U.S. Pat. No. 6,687,659 (incorporated by reference in its entirety herein) teaches a specific ABC of a certain class of ABCs, known as “perfectly matched layers” (PML).
Another area in which ABCs are required is that of electromagnetic radiation, which is governed by Maxwell's equations. In fact, the PML method was originally devised for such problems; see [3,4]. A typical problem in this area consists of one or more sources of electromagnetic radiation and some objects reflecting the electromagnetic radiation. One then wishes to find the values of the electromagnetic radiation field at grid points in some volume of interest. Since electromagnetic radiation extends to a potentially infinite domain, it is necessary to impose ABCs at the boundaries of the domain of interest.
In all of the above-mentioned areas, there are one or more field values to be determined in the volume of interest by solving the numerical problem. For example, in the Helmholtz equation, it is necessary to determine the value of the displacement, usually denoted by u, at each grid point. In electromagnetic radiation, there are two values to determine: the electric field and the magnetic field. PML methods specify certain equations that must hold in a layer of several extra grid points surrounding the original grid. Such a layer increases the size of the grid very significantly, and consequently, the number of equations is also increased. In contrast, the ABCs that use directional derivatives require just one extra layer of grid points, and thus they are more efficient than PML methods in terms of the computer time and memory required to solve a given problem.
The foundation for several ABCs is a certain mathematical condition, known as the Sommerfeld radiation condition [12]. This radiation condition involves the directional derivatives of the field value u, with the direction being away from the source of impact causing the waves. This direction is called the radial direction. The radiation condition provides a mathematical limit condition which must hold for the Helmholtz equation. ABCs that are based on the Sommerfeld condition strive to mimic the Sommerfeld radiation condition on the boundary of the grid.
Two well-known examples of ABCs that are based on the Sommerfeld radiation condition are the Engquist-Majda (EM) ABCs [5] and the Bayliss-Gunzburger-Turkel (BGT) ABCs [2]. Both of these methods have certain limitations. The EM ABCs use directional derivatives of the field value u in the direction that is orthogonal to the boundary of the domain. On box-shaped domains (which are widely used in practice), the orthogonal direction can deviate by almost 90° from the radial direction. Additionally, the orthogonal direction is not continuous across the edges and corners of such box-shaped domains. It is well known in the literature that the continuity problem contributes very significantly to inaccuracies in the solution of the Helmholtz equation, and there have been some attempts to handle this problem, such as [1]. However, none of the proposed methods offer the simple solutions presented by the Gradient Method.
The BGT method, as it was initially disclosed in [2] had certain configuration restrictions: the domain was assumed to be a sphere, and the source of impact was assumed to be at the center. Reference [9] presents an extension of BGT to elliptical bodies; however, as noted above, box-shaped domains are used extensively in many situations.
U.S. Pat. No. 6,687,659 to Shen discloses a hybrid approach for implementing absorbing boundary conditions in three dimensional (3D) finite difference (FD) acoustic applications such as post-stack and pre-stack seismic migration, and forward modeling.